library(tidyverse)
library(readr)
library(ProbBayes)



buffalo <- read_csv("data/buffalo_snowfall.csv")


data <- buffalo[59:78,c('SEASON', 'JAN')]

ybar <- mean(data$JAN)
print(nrow(data))
print(sd(data$JAN))
se <- sd(data$JAN) / sqrt(nrow(data))

post1 <- normal_update(c(10, 3), c(ybar, se))

many_normal_plots(list(c(10,3),
                       c(ybar,se),
                       c(post1[1], post1[2]))) +
  theme(legend.position = "none") +
  theme(text=element_text(size=18)) +
  xlab(expression(mu)) +
  annotate(geom="text",x=10,y=0.15, label="Prior", color="blue", size=7) +
  annotate(geom="text",x=ybar,y=0.10, label="Likelihood", color="red", size=7) +
  annotate(geom="text",x=post1[1],y=0.15, label="Posterior", color="pink", size=7) +
  ylab("Density")


x <- seq(-10, 10, length = 500)
y <- dcauchy(x, location = 0, scale = 1)
plot(x, y, type = "l", main = "Standard Cauchy Distribution")


x <- seq(-10, 10, length = 500)
mean <- 0
sd_list <- c(0.5, 1, 5)

plot(x, dnorm(x, mean, sd_list[1]), type = "l", col = "red", lwd = 2,
     ylab = "Density", main = "Normal Distributions with Same Mean")
for (sd in sd_list[-1]) {
  lines(x, dnorm(x, mean, sd), col = rainbow(length(sd_list))[which(sd_list == sd)], lwd = 2)
}
legend("topright", legend = paste("SD =", sd_list),
       col = c("red", rainbow(length(sd_list) - 1)), lwd = 2)













